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We calibrate an effective-one-body (EOB) model to numerical-relativity simulations of mass ratios 
1,2,3,4, and 6, by maximizing phase and amphtude agreement of the leading (2,2) mode and of 
the subleading modes (2,1), (3,3), (4,4) and (5,5). Aligning the calibrated EOB waveforms and 
the numerical waveforms at low frequency, the phase difference of the (2, 2) mode between model 
and numerical simulation remains below ~ 0.1 rad throughout the evolution for all mass ratios 
considered. The fractional amplitude difference at peak amplitude of the (2, 2) mode is 2% and 
grows to 12% during the ringdown. Using the Advanced LIGO noise curve we study the effectualness 
and measurement accuracy of the EOB model, and stress the relevance of modeling the higher-order 
modes for parameter estimation. We find that the effectualness, measured by the mismatch, between 
the EOB and numerical-relativity polarizations which include only the (2, 2) mode is smaller than 
0.2% for binaries with total mass 2O-2OOA/0 and mass ratios 1,2,3,4, and 6. When numerical- 
relativity polarizations contain the strongest seven modes, and stellar-mass black holes with masses 
less than 5OM0 are considered, the mismatch for mass ratio 6 (1) can be as high as 7% (0.2%) when 
only the EOB (2, 2) mode is included, and an upper bound of the mismatch is 0.5% (0.07%) when all 
the four subleading EOB modes calibrated in this paper are taken into account. For binaries with 
intermediate-mass black holes with masses greater than 5OM0 the mismatches are larger. We also 
determine for which signal-to-noise ratios the EOB model developed here can be used to measure 
binary parameters with systematic biases smaller than statistical errors due to detector noise. 



PACS numbers: 04.25. D-, 04.25.dg, 04.25.Nx, 04.30.-w 

I. INTRODUCTION 

Binary systems composed of black holes and/or neu- 
tron stars, spiraling in toward each other and losing en- 
ergy through the emission of gravitational waves, are 
among the most promising detectable sources of gravita- 
tional waves with the Laser Interferometer Gravitational- 
wave Observatory (LIGO) [1], Virgo [2], GEO [3], the 
Large Cryogenic Gravitational Telescope (LCGT) [4], 
and future space-based detectors. The detectors' noise 
level and the weakness of the waves prevent observing the 
waveforms directly. For this reason the search for gravi- 
tational waves from binary systems and the extraction of 
parameters, such as the masses and spins, are based on 
the matched-filtering technique, which requires accurate 
knowledge of the waveform of the incoming signal. 

The post-Newtonian (PN) expansion is the most 
powerful approximation scheme in analytical relativ- 
ity capable of describing the two-body dynamics and 
gravitational-wave emission of inspiraling compact bi- 
nary systems [5-8]. The PN approach expands the Ein- 
stein equations in the ratio of the characteristic velocity 
of the binary v to the speed of light or the characteristic 
size of the compact body to the relative distance between 
the two bodies. However, as the bodies approach each 



other towards merger, we expect the PN expansion to lose 
accuracy because the velocity of the bodies approaches 
the speed of light, and the relative distance becomes com- 
parable to the size of the compact body. The difficulty in 
analytically solving the Einstein equations in the merger 
regime lies mainly in its nonlinear structure. Solving the 
Einstein equations numerically overcomes this problem. 

Prior to the numerical-relativity breakthroughs [9-11], 
a new and unique method was proposed in analytical 
relativity to describe the dynamics and gravitational- 
wave emission of binary black holes during inspiral, 
merger and ringdown: the cffectivc-onc-body (EOB) ap- 
proach [12-16]. This approach uses the very accurate 
results of PN theory. However, it does not use those 
results in their original Taylor-expanded form (i.e., as 
polynomials in v/c), but instead in some appropriate re- 
summed form. In particular, the effective-one-body ap- 
proach [12, 14, 15, 17, 18] maps the dynamics of two 
compact objects of masses mi and m2, and spins Si 
and S2, into the dynamics of one test particle of mass 
/i = mim2/(mi -|- 7712) and spin moving in a de- 
formed Kerr metric with mass M = mi + m2 and spin 
Skoit- The deformation parameter is the symmetric mass 
ratio mim2/(mi -|- 7712)^ which ranges between (test 
particle limit) and 1/4 (equal-mass limit). The other 
crucial aspect of the EOB approach is the way it builds 
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the full waveform, including merger and ringdown. The 
EOB approach assumes that the merger is very short in 
time, although broad in frequency, and builds the merger- 
ringdown signal by attaching to the plunge signal a su- 
perposition of quasinormal modes. This match happens 
at the EOB light ring (or photon orbit) where the peak of 
the potential barrier around the merged black hole sits. 

The analyses and theoretical progress made in 
Refs. [18-32] have demonstrated that it is possible to 
devise and calibrate analytical EOB waveforms for use 
in detection searches. This is crucial, since thousands of 
waveform templates need to be computed to extract the 
signal from the noise, an impossible demand for numer- 
ical relativity alone. For example the EOB waveforms 
calibrated to numerical-relativity waveforms in Ref. [20] 
have been used in LIGO and Virgo to search for the first 
time for high-mass merging black holes [33]. 

This paper is a step forward in building more faithful 
EOB waveforms to be used for detection and parame- 
ter estimation. We calibrate the EOB model to accurate 
numerical-relativity simulations of mass ratios 1,2,3,4, 
and 6, so that the phase and amplitude agreement of the 
leading (2, 2) mode and also the subleading modes (2, 1), 

(3.3) , (4,4) and (5,5) are minimized throughout inspi- 
ral, merger and ringdown. The numerical simulations 
are produced by the pseudospectral code SpEC of the 
Caltech-Cornell-CITA collaboration [34-40] (see particu- 
larly Ref. [41] for details). The waveforms are extracted 
as Regge- Wheeler- Zerilli data [42], and extrapolated to 
infinite extraction radius [43]. Since the numerical- 
relativity modes satisfy the relation him = (— l)^/i|_m 
with high accuracy, where * denotes complex conjugate, 
we assume its validity also for the analytical modes. As a 
consequence, any statement in the paper concerning an 
(£, m) mode automatically holds for its complex conju- 
gate {£, —m) mode. We find that the (3, 2) mode has a 
distinct feature that currently cannot be accounted for 
with the EOB model used in this paper. Moreover, we 
find that the (6, 6) mode (and very likely other modes 
with m ~ t) can be calibrated in the same way as the 

(4.4) and (5,5) modes. However, we do not consider the 
(6, 6) mode since, for the range of mass ratios consid- 
ered here, its amplitude is much lower than the other 
subleading-mode amplitude, and therefore contributes 
little to the full polarization waveforms. 

The paper is organized as follows. In Sec. II, we de- 
scribe the EOB dynamics, the waveforms and its ad- 
justable parameters. In Sec. Ill we discuss the numerical- 
relativity simulations produced by the pseudospectral 
code SpEC [41] and estimate the phase and amplitude 
errors. Then, we calibrate EOB to numerical-relativity 
modes and discuss its effectualness and measurement 
accuracy when searching for gravitational waves with 
Advanced LIGO detectors. In Sec. V we compare our 
EOB model and its performance in matching numerical- 
relativity results to previous work. We summarize our 
main conclusions in Sec. VI. In Appendix A we discuss 
some interesting features of the (3, 2) mode. In Appendix 



B we list several quantities which enter the EOB wave- 
forms and energy flux. 



II. EFFECTIVE-ONE-BODY MODEL 

In Sees. IIA-IIC (see also Appendix B) we shall dis- 
cuss in detail all the building blocks of the EOB dy- 
namics and waveforms, and its adjustable parameters. 
The EOB model used in this paper is presented in a 
self-contained way to allow readers to reproduce it if 
desired. We note that many important features of the 
EOB model have been developed in several papers [12- 
14, 19, 20, 22, 23, 27-32, 44, 45]. 

If the EOB model were compared to the numerical- 
relativity simulations used in this paper without any cal- 
ibration, i.e., at the PN order currently known, 3PN in 
the conservative dynamics and 3.5PN in the radiation- 
reaction sector, we would find at merger a phase differ- 
ence for the (2,2) mode of up to 3.6 rad over the mass- 
ratio range q = 1,2,3,4,6. Moreover, the EOB am- 
plitude would peak around 30M before the numerical- 
relativity peak, with a fractional amplitude difference at 
the peak of ^ 8%. A straightforward way of reducing 
the differences is to insert in the dynamics and radiation- 
reaction force higher-order (pseudo) PN terms (or EOB 
adjustable parameters) and calibrate them to the numer- 
ical results. The advantage of the EOB approach is that 
the dynamics and radiation-reaction force (and modes) 
arc written in a way which isolates the crucial functions 
that determine the evolution. As we shall see below, 
these functions are the EOB radial potential A{r), or 
time-time component of the EOB metric, and some phase 
and amplitude functions appearing in the EOB (factor- 
ized) gravitational modes. 



A. EfFective-one-body dynamics 

We set M = mi -I- 7712, = mim2/M = vM, q = 
mi/m2, and use natural units G = c = 1. The EOB 
effective metric reads [13] 

dsls = -~A{r) dt^ + ^ dr^ + r^ (dO^ + sin^ 6 , 

(1) 

where we use dimensionless polar coordinates (r, $) and 
their conjugate momenta {pr,p<^). Replacing the radial 
momentum pj. with p^^ which is the conjugate momen- 
tum to the EOB tortoise radial coordinate r* , 

dr, ^ /D(0 

dr A{r) ' ^ ' 
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we obtain the EOB effective Hamiltonian [13, 14, 



(3) 



where we have neglected the factor Z?(r)^/A(r)"' in front 
of the term p^^ which would introduce PN terms higher 
than 3PN order, but more importantly would cause the 
EOB gravitational frequency to grow too quickly near 
merger. 

The real EOB Hamiltonian reads [13] 



= M ^ l + 2v 



- M . (4) 



The Taylor approximants to the coefficients A{r) and 
D{r) can be written as [13, 14] 



(5) 



i=0 



i=0 



The functions A{r), D{r), Ak{r) and Dk{r) all depend 
on the symmetric mass ratio v through the z^-dependent 
coefficients aj;(z^) and di{i') [see Eqs. (47) and (48) in 
Ref. [22]]. The functions Ak{r) and Dk{r) are currently 
known through 3PN order, i.e., fc = 3. During the last 
stages of inspiral and plunge, the EOB dynamics can be 
adjusted closer to the numerical simulations by includ- 
ing in the radial potential A(r) a few adjustable parame- 
ters of the EOB dynamics. Notably, the 4PN coefficient 
a5(i/) [20, 22, 23, 27-29, 44] and even the 5PN coefficient 
aeii^) [31]. 1 

To enforce the presence of the EOB innermost stable 
circular orbit (ISCO), Ref. [14] suggested using the Pade 
expansion of the function A{r). For A{r) we employ the 
Pade expression Al{r) at 5PN order, while for D{r) we 
use the Pade expression D^{r) at 3PN order. We could 
also introduce EOB adjustable parameters at 4PN and 
5PN order in D{r), say ^4(1^) and ^5(1^). However, this 
modification would affect mainly the radial motion [see 
Eq. (10a) below] which becomes important only at the 
very end of the evolution. For the EOB model devel- 
oped in this paper we find that these other adjustable 
parameters arc not needed. The quantity -D§(r) reads 



while A'!^{r) reads 



(52i^-6i/2) + 6i/r + r3 



Num(A^) 
Den(^i) 



(6) 



(7) 



EOB dynamics 


EOB waveform 


adjustable parameters 


adjustable parameters 


as, ae 


"''match 












pQNM 



TABLE I: Summary of adjustable parameters of the EOB 
model considered in this paper. We notice that to calibrate 
the EOB (2,2) mode, we only need asjOg and At^^^^^^. To 
calibrate each subleading mode (2, 1), (3, 3), (4, 4), and (5, 5), 
we need four adjustable parameters. The values of the ad- 
justable parameters used in this paper are given in Eqs. (36) 
to (39) and (41). 



with 

Num(.4^) 
and 



= r'* [-64+ 12a4 4a5 + ae 
+ [32 - 4 04 — 05 — 24i/] , 



64i/ - Aiy^ 



(8) 



32 t/^ 



(9) 



Deii{Al) = 4 04 + 4 04 05 + flg — 04 ae + 16 ag + (32 04 

-f 16 as - 8 ae) + 4 04 + 32 + r [4 04 4- 04 05 

-|-16 + 8 ag + (32 04 — 2 ae) + 32 + 8 

[16 04 + 8 05 + 4 ag + (8 a4 + 2 05) v 

[8 04 4 ag + 2 ag + 32 1/ - 8 v"^] 

+r* [4 04 + 2 05 + ag + 16 1/ - 4 1/^] 

[32 - 4 a4 - as - 24 i/] , 

where = [94/3 — (41/32) tt^] 1/ and to ease the nota- 
tion we have omitted the v dependence of as and ag in 
the expressions above. The quantities as and ag are the 
adjustable parameters of the EOB dynamics [23] (see Ta- 
ble I). They will be determined below when calibrating 
the EOB to numerical-relativity waveforms. Their ex- 
plicit expressions are given in Eq. (36). 

The EOB Hamilton equations are written in terms of 
the reduced, i.e., dimensionless quantities iJ™*^' [defined 
in Eq. (4)] [12]. They read^ 



dr 
dt 
d$ 
dt 

dpr, 

dt 

dp^ 
dt 



Al{r) dH 



real 



dpr, 
{r,Pr,,P'S-) 



{r,Pr,,P^) 



(10a) 



dp<s, 
Al{r) 



(10b) 



nK 7- 



9?- 



ir,Pr,,P^>) 



nK 



Pr, 



(lOc) 

(lOd) 



^ We notice that the second term on the right-hand side of 
^ The radial potential A{r) may contain logarithmic terms at 4PN Eq. (10c) is generated when taking the nonspinning limit of the 

and 5PN orders [47, 48] which we do not try to model here. spinning EOB model of Ref. [16] 
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with the definition fl = d^jdt = M^. The initial 
conditions for tlie EOB Hamilton equations will be dis- 
cussed in Sec. II D. Furthermore, for the $ component of 
the radiation-reaction force we use a non-Keplerian (nK) 
radiation-reaction force:"^ 



(11) 



where = ^^1'^^ and dE/dt is the gravitational- wave 
energy flux for quasi-circular orbits obtained by summing 
over the gravitational- wave modes {l,m). We use 



dE 

m 



7 t 

E E 

1=2 rn=£-2 



(12) 



Note that because |/i^,m| = we extend the sum 

over positive m modes only. Moreover, for the cases stud- 
ied in this paper, including more modes in the summation 
has a negligible effect on the energy flux. Specifically, if 
we sum £ through ^ = 8 and sum m — 0,...,£, the 
gravitational- wave phase of the EOB (2, 2) mode changes 
by less than 0.01 rad at merger, which is negligible com- 
pared to the phase error of numerical-relativity modes. 
We find that the dominant computational cost in gener- 
ating the EOB waveforms is the calculation of the energy 
flux. The choice of modes {i,m) in Eq. (12) saves us 
about a third of the computational time when compared 
to the case where the sum extends up to ^ = 8, and runs 
over m = 0, . . . ,£. The explicit expression of the modes 
hem is given below, in Sees. II B and II C. 



B. EOB waveform: Inspiral & Plunge 

Having described the inspiral and plunge dynamics, we 
now turn to the gravitational- wave modes him- The lat- 
ter can be employed to compute consistently the inspiral 
dynamics through the radiation- reaction force [31] [see 
Eq. (12)]. The inspiral and plunge EOB modes are given 
by 



insp— plunge 



(13) 



where the Ni,n describe effects that go beyond the 
quasi-circular assumption and will be defined below [see 
Eq. (22)], and the hj^ are the factorized resummed 
modes. In the nonspinning case, Refs. [23, 31] have 
shown that the resummed, factorized modes proposed in 
Rcf. [30] are in excellent agreement with the numerical 
waveforms [29, 34]. We have [30], 



'''em '-'nff ^ «m e 



■'cff 



(a 



•em) 



(14) 



^ Note that Eq. (11) is only implicitly non-Kcplcrian, in contrast 
to similar expressions in other papers which explicitly introduce 
non-Keplerian terms. In this case, the non-Kcplcrian behavior is 
hidden in the wave amplitudes /t£„ , as described in the following 
section. 



where e denotes the parity of the multipolar waveform. 
In the circular-orbit case, e is the parity of £ -|- to: 



, £ -f TO is even 
1 , ^' + 771 is odd 



(15) 



The leading term in Eq. (14), h^^^ 
contribution 



is the Newtonian 



Mv 



n',ice+.{iy)ViY'-^'-^ (|,$) , (16) 



where TZ is the distance from the source; the F^™(0,$) 
are the scalar spherical harmonics; the functions and 
ce-\-t{v) are explicitly given in Appendix B [sec Eqs. (B7) 
and (B8)]. Moreover, for reasons that will be explained 
in Sec. HIE, wc choose 



U+e-2) 



(£,TO)^(2,1),(4,4), (17a) 
(€,to) = (2,1), (4, 4); (17b) 



with u$ and defined by [22] 

v^ = nrn = f)r [■(/'(r,p$)]^/^ 



(18) 



and 



7/'(r,_p$) 



[l + 2u yAl{r){l+pi/r^)-l]] 



dAl{r)/dr 



(19) 

The quantity p$ in the above equation is the dynamical 
p<i> being used in the evolution, contrary to the choice 
made in Ref. [22] where p$ was chosen to satisfy the 
circular-orbit condition. The function S'^g in Eq. (14) is 
an effective source term that in the circular-motion limit 
contains a pole at the EOB light ring. It is given in terms 
of the EOB dynamics as 



e = 0, 
e = 1, 



(20) 



where H'^^ {r,pr,,p^) can be read from Eq. (3). The fac- 
tor Tgm in Eq. (14) resums the leading order logarithms 
of tail effects, it reads 

Tim 7-„ ^ exp \TTmn 

T{£+1) J (21) 

X exp [2i m n i/""^' log(2 m n tq)] , 

where Tq = 2M/y/e [26] and can be read from 

Eq. (4). 

The factor e'-^'"^ in Eq. (14) is a phase correction due 
to subleading order logarithms, while the factor {pemY 
in Eq. (14) collects the remaining PN terms. The full 
expressions for dem and pem are given in Appendix B. 
To improve the agreement with the numerical-relativity 
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FIG. 1: Amplitude of extrapolated numerical-relativity wave- 
forms for the dominant modes. The two panels from top to 
bottom are for mass ratios q — 1 and q = 6, respectively. 
Each curve is labeled with its respective {I, m) mode, and the 
time-delay /^t^^y^ between the extrema of |/i22| and |fe!m|. 
The horizontal axis measures the time-difference to the peak 
of 1/122]. 



higher modes, we introduce and cahbrate in 5im and 
Pirn a few higher-order, yet unknown, PN terms (see 
also Refs. [49, 50]). Details of this are given below [see 
Eqs. (38) and (39)], so here we merely remark that in 
the spirit of Ref. [23] , those coefficients should be consid- 
ered adjustable parameters of the EOB waveform (sec Ta- 
bic I). The introduction of these higher-order PN terms in 
the modes ^ (2, 2) is not surprising, because these modes 
are known at lower PN order than the (2,2) mode. We 
note that the adjustable parameters in 5im and pim will 
be used only to improve the agreement between the EOB 
and numerical-relativity modes. They will not be in- 
cluded in the energy flux entering the dynamics through 
Eq. (12). 



Finally, the function entering Eq. (13) reads 



Pr' 



p3/2 



X cxp 



Pr' 



Pr' 

rCl 



(22) 



where the quantities a- and 6^ are the non- 
quasicircular (NQC) orbit coefficients [23, 25, 27-29, 32, 
45]. 

To better understand how we fix the parameters 
and 6^"", we plot in Fig. 1 the amplitudes of 
the dominant numerical him- As already observed in 
Refs. [20, 32, 51-53], the peaks of the modes occur at 
different times. In Fig. 1 we have indicated these times 
relative to the peak of /i22- Our goal is to model the EOB 
modes, through the parameters a^''^'" and 6^'^'" , in such a 
way to (i) reproduce the shape of the numerical-relativity 
amplitudes close to the peak, and (ii) preserve the time 
differences between the modes. 

So, for each mode, we fix the three parameters a^*"' in 
the amplitude and the two coefficients fc^**™ in the phase 
by requiring that the peaks of the numerical and EOB 
as well as their frequencies at the peaks, coincide.^ 
Specifically, we impose five conditions listed below for 
each mode. 

1. The time at which the EOB ^,22 reaches its peak 
should coincide with the time at which the EOB or- 
bital frequency reaches its peak. We denote this 
time with t^^^j^- It was observed [23, 26, 31] that, 
once the EOB and numerical phases are aligned at 
low frequency and calibrated, the time at which 
the numerical /i22 reaches its peak coincides with 
the EOB light-ring time. Moreover, the latter oc- 
curs immediately (< IM in time) after the peak 
of fl. The peaks of higher-order numerical modes 
differ from the peak of the numerical /122 mode by 
a few M in time. We define this time difference as 



At 



peak 



peak peak 



22 



— ^d\h^^\/dt=0 - td\h^^\/dt=0 1 

and require that the peaks of the EOB him occur 



(23) 



at the time t^^^i. 



2. The peak of the EOB him should have the same 
amplitude as the peak of the numerical him, that 
is 



I'^m ('■peak + ^"-peakJ I ~ | ''f m C-poak J | 



(24) 



During the course of this work we noticed that Rof. [32] had inde- 
pendently developed the same procedure of us in fixing the NQC 
phase coefficients using io^^' and However, Ref. [32] com- 

putes those numerical quantities at the light-ring position 3M 
instead of the peak of the (£, m) modes as we do. Reference [32] 
focused on the EOB modeling in the extreme mass-ratio limit. 
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FIG. 2: We compare the numerical-relativity and BOB /i22 amplitudes with and without the NQC corrections Nirn given in 
Eq. (22). We also plot the numerical and EOB gravitational frequency of the (2, 2) mode and twice the EOB orbital frequency. 
The left panel refers to g = 1 and the right panel to g = 6. The horizontal axis is the retarded time in the numerical-relativity 
simulation. The vertical lines mark the peaks of the numerical-relativity /122 amplitudes. 



3. The peak of the EOB him should have the same 
second-order time derivative as the peak of the nu- 
merical h(rm that is 



(25) 



This condition guarantees that the local extremum 
of I hfo^ \att = t^^^^ + At^Zl^ is a local maximum. 

4. The frequency of the numerical and EOB him wave- 
forms should coincide at their peaks, that is 



^°^(*poak + ^4"ak) — '^^^(^pTak) (26) 



5. Time derivative of the frequency of the numerical 
and EOB him waveforms should coincide at their 
peaks, that is 



UJi, 



EOB/,0 

l''pcak 



peak 



'-'em V^pcak/ 



(27) 



In Sec. 



III A 

peak'' 



shall find 
d'\hZ\/dt'\ 



that the functions 



,NR 



(-lira \ 



and 



(^poak) ^'"'5 reasonably approximated by smooth func- 
tions of V (see Table III). 

The NQC coefficients a^'" are calculated within the 
EOB model using the fitting formulas in Table III. The 
calculation involves a computationally expensive itera- 
tive procedure. Basically, in each round of the iteration, 
g^fm g^jjj ^'ifm j^j^g calculated to satisfy the five conditions 

listed above. The amplitude corrections a^*"" then enter 
the dynamics through the energy flux given in Eq. (12). 
The new EOB dynamics generate new modes and, thus, 
new a^'"*. We stop the iteration when a^'"" converge 
on successive runs. In this paper, for technical conve- 
nience, we include only al^^{v) in the energy flux (12) 



and ignore the effect of higher-order-mode NQC correc- 
tions on the inspiral dynamics. However, we do include 
the higher-order-mode NQC corrections when building 
the EOB waveforms and compare them to the numerical- 
relativity ones. Neglecting higher-order-mode NQC cor- 
rections in the energy flux (12) is insignificant for three 
reasons: The higher-order-mode contribution to the en- 
ergy fiux is about an order of magnitude smaller than 
that of the dominant (2, 2) mode; the NQC corrections 
in the amplitude are a relatively small correction, typi- 
cally ^ 10% at merger; and the NQC correction is most 
important close to merger where the radiation reaction 
has little effect on the plunging dynamics. 



The iterative procedure explained above usually takes 
4 to 5 iterations to converge, bringing a factor of 4 to 5 
to the computational cost of generating an EOB wave- 
form. Therefore, to reduce computational cost, we give 
fitting formulas of al^^ in Eq. (40). We shall emphasize 
however that a"== are not adjustable parameters of the 
EOB model and are not required as inputs to generate 
the (2, 2) modes. 



Finally, to demonstrate the effect of the NQC correc- 
tions (22), we show in Fig. 2 the EOB amplitude for the 
(2,2) mode without and with NQC corrections, that is 
without and with the factor Nim in Eq. (13). We notice 
that for q = 1 and g = 6, the amplitude without the 
NQC corrections differs from the numerical one only by 
-0.23% and -0.67%, at t = 6.2Af and t = 5.7M before 
the peak, respectively. However, even this small differ- 
ence needs to be removed to minimize the error when at- 
taching the merger-ringdown waveform. The latter will 
be discussed in the next section. 
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C. EOB waveform: Merger & Ringdown 

The procedure of building the merger-ringdown wave- 
form in the EOB approach has improved over the 
years [12, 19, 20, 23, 28, 29, 54]. For each mode (£, m) 
we have 



N-l 



^(t-*L"tch) 



(28) 



where n is the overtone number of the Kerr quasinormal 
mode (QNM), N is the number of overtones inchided 
in our model, and Aimn arc complex amplitudes to be 
determined by a matching procedure described below. 
The quantity ct^,„„ = uiimn — i/Timn, where the oscillation 
frequencies iVimn > and the decay times Ti„in > 0, 
are numbers associated with each QNM. The complex 
frequencies are known functions of the final black-hole 
mass and spin and can be found in Ref. [55]. 

In this paper we model the ringdown modes as a lin- 
ear combination of eight QNMs, i.e., N ~ 8. Mass and 
spin of the final black hole Mf and a f are computed from 
numerical data for mass ratios q = 1,2,3,4, and 6. No- 
tably, we employ the fitting formula obtained by fitting 
the numerical results of Mf and af 



M 




0.4333i^2 - 0.4392^^^(29a) 



^ = \/T2iy - 3.8711/2 + 4.028i^^. 



(29b) 



The above formula differs from the analogous fitting for- 
mula given in Ref. [20] by < 0.3% in Mf and < 2% in 
a f , because of the more accurate numerical data used in 
this paper. 

The complex amplitudes Aimn in Eq. (28) arc deter- 
mined by matching the EOB merger-ringdown waveform 
(28) with the EOB inspiral-plunge waveform (13). In 
order to do this, N independent complex equations are 
needed. In Ref. [23] we introduced the hybrid-comb 
matching in which TV equations are obtained at — 4 
points evenly sampled in a small time interval At^^^^ 
ended at i^tch' we imposed the condition that the 
inspiral-plunge and merger-ringdown waveforms coincide 
at the TV — 4 points and their first and second order time 
derivatives coincide at the first and the last points. Un- 
like in Ref. [23], we now no longer require second or- 
der time derivatives of the waveforms to coincide any- 
where in order to improve the numerical stability of the 
matching procedure. Instead, we impose the continu- 
ity of the waveform at TV — 2 points evenly sampled 
to t^^tch> and we require conti- 



irom Eniatch 

nuity of the first time derivative of the waveforms at 

j-fm _ \4-tm „ J -i-lm ■ „ 
''match "''match "^^^^ ''match i ^•'^•J 



1 insp- plunge /,£m 
'Vm v'-match 



k 



TV- 3 



"''match/ 



7 merger - 



-RD 



match 



TV- 3 
(fc = 0,l,2,.. 



-At 



match/ 



,TV-3), (30) 



and 



; insp— plunge / ,gm 
"'(.m I'match 



TV - 3 



"''match/ 



' merger— RD / 



V match 



TV-3 



match/ ' 



(fc = and TV - 3) . (31) 



The matching time tf"atch fixed to be the peak of the 
EOB h,„^ mode, i.e., 4",tch - ipeak + At^p^ak- The match- 
ing interval Ai^j.^[^ is an adjustable parameter that we 
fix by reducing the difference against numerical merger- 
ringdown modes (see Table I). 

Finally, the full (inspiral-plunge-mcrger-ringdown) 
EOB waveform reads 



''e{t 



match 



, mergor-RD n(, ,im 



match/ ■ 

(32) 

It was noticed in Ref. [23] that when the lowest QNM 
frequency is substantially larger than the EOB mode fre- 
quency at i^tch' the EOB mode frequency will generally 
grow very rapidly to the QNM frequency immediately af- 
ter t^tch- Such growth in the EOB frequency is much 
more rapid than what is seen in numerical-relativity fre- 
quencies. Moreover, when this happens, the EOB am- 
plitude shows an unphysical "second peak" shape where 
the ringdown amplitude grows for a while before eventu- 
ally decaying. The growth of the EOB frequency can be 
slowed down by including a pseudo QNM [23] that has 
a frequency close to the EOB mode frequency at i^t^h 
and a decay time comparable but smaller than the decay 
time of the least damped n ~Q QNM. As we shall discuss 
below, we find it necessary to introduce a pseudo QNM 
in modeling the EOB (4, 4) and (5, 5) modes. The pseudo 
QNM should be counted as another adjustable parame- 
ter of the EOB waveform (see Table I). The frequency 
and decay time of this pseudo QNM mode are given in 
Eq. (41). 

We have outlined the procedure to match the inspi- 
ral waveform to the merger-ringdown waveform. We 
would like now to understand what is the intrinsic error 
that this procedure introduces. To answer this question, 
we build an inspiral-merger-ringdown waveform where 
the inspiral part coincides with the numerical-relativity 
waveform, and the merger-ringdown part is built using 
the EOB procedure. We then extract the intrinsic error 
by comparing it to the numerical-relativity full waveform. 
The results for the (2,2) and (4,4) modes are shown in 
Fig. 3. By construction, the two waveforms agree ex- 
actly before the matching time i^tch' i-^-' the time of 
the peak amplitude. For /122, the relative amplitude dif- 
ference and phase difference during ringdown are about 
10% and 0.1 rad. These are reasonable intrinsic errors 
for the EOB model and are comparable to systematic er- 
rors in the best existing analytical models [23, 31]. For 
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FIG. 3: Amplitude (in units of MjTV) and frequency (in units of comparison between the full "NR" waveform and the 

"NR+QNM" waveform generated by attaching QNMs to the inspiral-plunge numerical waveform. We show also the relative 
amplitude and phase differences. In the left panel, we compare /i22. In the right panel, we compare the numerical /i44 mode 
with two "NR+QNM" mode. One of them is generated by attaching the physical QNMs, the other is generated by attaching 
both the physical QNMs and the pseudo QNM. The former is very different from the numerical-relativity mode and we do not 
show their amplitude and phase differences. All /i44 amplitudes have been multiplied by a factor of 20, so that they are more 
visible. The horizontal axis is the retarded time in the numerical-relativity simulation. 



/144, the pseudo QNM reduces the amplitude and phase 
differences substantially to the level of 50% (when the 
ringdown amplitude is below 10% of the peak amplitude) 
and 0.6 rad. Although the differences are not as small 
as those of /122, they are for now acceptable considering 
the relatively small amplitude of ^144 compared to at 
least for the mass ratios considered in this paper. So, in 
the following we shall not attempt to over-calibrate the 
EOB /144 model to obtain smaller differences against the 
numerical results. 

The intrinsic error depends on the procedure to match 
the inspiral waveform to the merger-ringdown waveform 
In particular, it depends on the choice of t 
At 



match 



match' 



as well as the continuity conditions we impose 
on the points sampled from t^.^j^ - At^^^,^^ to t^™^^^. In 
Fig. 3, the results are optimized only over At'^^^^. To 
find the best matching procedure, in principle, we should 
optimize over both tf™tch ^^mStch' ^^"^ consider op- 
tions that sample points and impose continuity condi- 
tions in different ways. Since the relation between the 
intrinsic error and all these parameters is not straightfor- 
ward, we decide to refrain from fine-tuning the ringdown 
waveform by assigning different matching procedures to 
different mass ratios or modes. We prefer to use a single, 
simple prescription that works well for all mass ratios and 
modes. Thus, we fix tf™^^!^ to be the peak of the EOB 
him waveform and find that this prescription works fairly 
well. 



In summary, there is an intrinsic error introduced 
by the current procedure to match inspiral to merger- 







g = 2 


g = 3 


g = 4 


g = 6 


U/M 


820 


770 


570 


670 


870 


t2/M 


2250 


2255 


1985 


1985 


2310 



TABLE II: The range of integration (ti,t2) for waveform 
alignment and mass ratios q = 1, 2, 3, 4 and 6. 



ringdown EOB waveforms. This error cannot be im- 
proved by better calibrating the EOB inspiral-plunge dy- 
namics and waveforms. It can be overcome only by im- 
proving and/or changing the matching procedure. We 
leave to the future this important work. 



D. Initial conditions for the EOB dynamics 

Before completing this section, we briefly review the 
way in which initial conditions of the EOB Hamilton 
equations (10) are implemented. 

In Ref. [16], quasi-spherical initial conditions are given 
for generic precessing black hole binaries. We adopt the 
nonspinning limit of those conditions 



(33a) 



g^rcal 


= 0, 




dr 




^^rcal 


1 dE 


{d^H'^^^drdp^) 


dPr. 


V dt 


(5^rcal /Qp^ ) (aSiJrcal / Q^2 ) 


^^rcal 


= ^0- 









,(33b) 



(33c) 
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Given the initial orbital frequency flo, we solve the initial 
r, and p$ from Eqs. (33a)-(33c) and set the initial or- 
bital phase $ = 0. We always start the orbital evolution 
at initial orbital frequency MOg < 0.0025, correspond- 
ing to an initial radius r > 50M, such that the binary 
orbits are sufficiently circularized at the frequency where 
numerical simulations start. ^ 



III. CALIBRATIONS AND COMPARISONS 

In this section we calibrate the EOB model against the 
numerical simulations, and compare numerical and EOB 
(2,2), (2,1), (3,3), (4,4) and (5,5) modes. 



A. Numerical-relativity simulations of 
unequal-mass binary black holes 

The numerical simulations themselves are described in 
a separate paper [41]. We extract both the Newman- 
Penrose (NP) modes ^'4™ and the strain modes him from 
the simulations. The strain modes are extracted with 
the Regge-Wheeler-Zerihi (RWZ) formalism [42, 56-58] 
(see Appendix A of Ref. [23] for details of the numerical 
implementation used to obtain him)- The waveforms are 
then extrapolated to infinite extraction radius with order 
N = 5 polynomials in the q = 1 case, and TV = 3 poly- 
nomials in other cases [43]. In this section, we use the 
RWZ him to calibrate the EOB adjustable parameters 
and to determine the EOB NQC coefficients as functions 
of the mass ratio. We use the NP ^'4'" only to align nu- 
merical and/or EOB waveforms at low frequency where 
numerical errors of the former are smaller than those of 
the RWZ hem- 

We adopt the same waveform alignment procedure 
used in Refs. [22, 23, 26]. That is, we align waveforms at 
low frequency by minimizing the quantity 

E{At, A0) = / ' [Mt) -M-t- At) - A0]2 dt , (34) 

over a time shift At and a phase shift A^, where 0i(t) 
and 4>2{t) are the phases of the two ^£'4™ waveforms. The 
range of integration {ti,t2) is chosen to be as early as 
possible to maximize the length of the waveform, but 
late enough to avoid the contamination from junk radia- 
tion present in the numerical initial data. The range of 
integration should also be large enough to average over 
numerical noise. Moreover, this range should extend from 
peak to peak or trough to trough of oscillations (if visi- 
ble) in the gravitational-wave frequency due to residual 



^ Note that Eq. (33b) is derived in Ref. [16] for pr- At large initial 
separations, the difference between pr» and pr is negligible and 
we do not distinguish them when setting initial conditions. 
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Fit formula 
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10.67 - 41.41 -h 76.1 
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(5m (0.1867 + 0.6094 1/) 
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-vSm (0.09051 - 0.1604^) 
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2) 


0.2733 + 0.2316 v + 0.4463 u'' 






(3 


3) 


0.4539 + 0.5376 v + 1.042 v'^ 


^em l^p 


m \ 
cak/ 


(2 


1) 


0.2907 - 0.08338 v + 0.587 v'^ 






(4 


4) 


0.6435 - 0.05103 v + 2.216 v'^ 






(5 


5) 


0.8217 + 0.2346 u + 2.599 u'^ 






(2 


2) 


0.005862 + 0.01506 v + 0.02625 i/^ 






(3 


3) 


0.01074 + 0.0293 v + 0.02066 v'^ 


^Im I '-peak j 


(2 


1) 


0.00149 + 0.09197 v - 0.1909 i/^ 






(4 


4) 


0.01486 + 0.08529 v - 0.2174 j/^ 






(5 


5) 


0.01775 + 0.09801 v - 0.1686 v"^ 



TABLE III: We list in the third column //-fits of At^ak, 

j.NRi'.fm \\ j2 |> NRI / 7,2 1 NR^.fm x „„J,-,NR/.fm \ 

"'£m v^'pcak cak/ 5 "j^i^ ^''pcak/ 

peak 

against numerical data of mass ratios q = 1, 2, 3, 4 and 6 for 
modes (2, 2), (3, 3), (2, 1), (4, 4) and (5, 5). In the fitting for- 
mulas, the relative mass difference is Sm = (mi — 1712) /{mi + 
m2) = VI — 4 !/. 



eccentricity in the initial data. For different mass ratios, 
the lengths of numerical simulations and waveforms are 
different, thus we must also choose different ti and t2 in 
Eq. (34). Ref. [59] suggested a minimal length of the in- 
tegration interval (^1,^2), which is satisfied by our choices 
as listed in Table II. 

The numerical uncertainties are estimated by combin- 
ing convergence errors with extrapolation errors. The 
convergence estimates also use the matching procedure 
described above. Specifically, the resolution convergence 
uncertainty is obtained by matching data from the high- 
est and second-highest resolutions over the relevant re- 
gion from Table II, then taking the relative amplitude 
and phase difference between them. The same process 
is repeated for the extrapolation uncertainty, comparing 
waveforms extrapolated with two orders — the order used 
in this paper, and the next higher order. Specifically, 
those orders are N = 5 and N = 6 for the q = 1 case, 
and iV = 3 and = 4 for g = 2, 3, 4, 6. The absolute val- 
ues of those uncertainties are then added to give the final 
uncertainty estimate. These uncertainties are shown as 
dotted lines in every figure of this paper where phase and 
fractional amplitude differences are shown. 
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FIG. 4: We show the amplitude of the numerical-relativity 
/i22 for mass ratios 5=1,2, 3, 4, 6. We have time shifted the 
modes so that their peaks are aligned. We have also rescaled 
them by v. The horizontal axis is the retarded time in the 
numerical-relativity simulation. The inset shows an enlarge- 
ment of the merger region. 



crtheless, we find the fitted values accurate enough for 
constraining the shape of the amplitude peaks. In fact, 
the error in the EOB merger-ringdown waveform of these 
modes are dominated by the intrinsic error of the match- 
ing procedure discussed in Sec. II C. The error on Atp™^^]^ 
is the error in determining the peak of \hi„i\ in the EOB 
model. Since the peak of the orbital frequency is used as 
the reference time of merger in the EOB model, and since 
it coincides with the peak of the numerical amplitude /i22 
to within l.SAf (see Sec. II B), the 0.5M error in fitting 
^^poak ^^'■^ ^'•peak Sufficiently small. The relative er- 
rors in fitting all other quantities are within 1% and we 
expect these fitting formulas to work with such accuracy 
for any mass ratio g < 6. 

It is interesting to observe that the reason why the 
shape and characteristics of the numerical amplitude 
peaks can be easily reproduced by polynomials in v (see 
Table III) rests on the fact that when the modes for dif- 
ferent v are aligned at the peak and rescale by ;/, their 
peaks differ by less than 7% and the width at half peak 
by less than 18%. This is illustrated in Fig. 4, and was 
initially pointed out in Ref. [52]. 



C. Calibrating the EOB adjustable parameters 



B. Extracting information from 
numerical-relativity waveforms 

As discussed in Sees. II B. II C we need to extract 
specific information from the numerical data that we 
use to determine the NQC coefficients a^*^" and fef '"' in 
Eq. (22), and the QNM coeflacients At^n in Eq. (28). 

In Table III we have listed the fitting formu- 
las for the relevant quantities At^p™,, \h'^^{t%ly)l 
dP\hZ\ldt\^ , u:ZK^.^.). and a;NR(t^™ J. These 

peak 

formulas are least-square fits of numerical-relativity re- 
sults at mass ratios q = 1,2,3,4 and 6 and numerical 
results in the test-particle limit v — 0.001 generated by 
a time-domain Teukolsky code [60]. The errors of the 
fitting formulas are worst for the (5, 5) mode and for 
the quantity (P\h^^\/dt^\ . The (5, 5) mode generally 

has the lowest amplitude among the modes being studied 
and is therefore mostly contaminated by numerical arti- 
facts. Specifically, we find oscillations in the amplitude 
of numerical waveforms that are amplified in the process 
of extrapolating the extraction radii to infinity. Such 
oscillations modify the position and shape of the peak 
amplitude and they start to become significant for the 
(5, 5) mode and the other modes with lower amplitudes. 
The errors suggest that we can barely model the merger 
amplitude of the (5, 5) mode for generic mass ratios. Nu- 
merical errors prevent us from modeling the merger and 
ringdown of any other mode with smaller amplitude. The 
relative error on <P'\h^^\/ dt^\^f^ may reach a few per- 
cent for the (2,2), (2,1), (3,3) and (4,4) modes. Nev- 



We carry out the calibration of the adjustable param- 
eters as follows. First, we fix the EOB-dynamics param- 
eters by minimizing the phase difference between the nu- 
merical and EOB (2, 2) modes during the inspiral. Sec- 
ond, we shall evolve the EOB dynamics with the cali- 
brated parameters and fixed the EOB-waveform param- 
eters by minimizing the difference between the numerical 
and EOB full waveforms of all relevant modes. Moreover, 
the adjustable parameters might be functions of the mass 
ratio v. So, we first calibrate them for individual mass 
ratios then fit the calibrated values with quadratic func- 
tions of v. We find that many parameters depend weakly 
on V and can be set as constants. 

Figure 5 summarizes the calibration result of the inspi- 
ral dynamics and our choice of the az{v) and 06(1^) values. 
We calibrate the EOB model to five numerical /122 wave- 
forms of mass ratios q = 1,2,3,4 and 6. For each mass 
ratio, we show a contour in the a^iv) / v-aQ{i') / v param- 
eter space in which the numerical and EOB /122 wave- 
forms agree in phase to within 0.2 rad at merger, that is 
at the peak of the (2, 2) mode. As observed in Ref. [31] 
(and also evident from Fig. 5), there is strong degeneracy 
between a^{v) and aQ{i'). Furthermore, there are many 
ways to choose 05(1/) and adiv) so that the numerical and 
EOB ft,22 waveforms agree for q~ 1, 2, 3, 4 and 6. For in- 
stance, we could choose a point near ac,[v)/v = — 10 and 
o.f){v) jv — 126 where the contours of different mass ratios 
approximately cross. However, we find that those values 
do not satisfy the constraint imposed by the self-force 
result of the ISCO shift of Ref. [61]. More importantly, 
even if we were not keen at imposing the ISCO shift pre- 
diction, the EOB model obtained by choosing ac^{v)/v 
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FIG. 5: We calibrate adjustable parameters of the EOB dy- 
namics. For mass ratios q = 1,2,3,4 and 6, the shaded re- 
gions correspond to (05,06) values for which the EOB and 
numerical-relativity /122 agree within 0.2 rad at merger, i.e., 
at the peak of the numerical h22- In the inserted subplot, for 
aQ{ii)/v = 184, we show ac,{u)/v values constrained by the 
shaded regions and by the test-particle ISCO-shift result [61], 
and also the quadratic fit (red curve) given by Eq. (36a). 



and a^(y)lv around that crossing region will not be very 
satisfactory, because the point will not lie in the middle 
of the contours for each mass ratio. 

By contrast, we follow a more satisfactory route. We 
decide to model as [v) / v and a^{v)lv s& generic quadratic 
functions of i/, i.e., ac,{v)/v = X]i=o "^i*^ ^'^'^ i6(i')/'^ = 



In the test-particle limit i/ — >■ 0, we impose 



the constraint on a'^^ and Og*''' that the conservative EOB 
dynamics incorporates the exact ISCO-frequency shift of 
the self- force calculation [61], that is [47] 



(35) 



and flg*^ are deter- 



The remaining five coefficients of 
mined by minimizing the distances of the a^lv) / v-ae^v)!/ 
points to the center of the contours in Fig. 5 (i.e., the 
points at equal distance from the top and bottom bound- 
aries of the contours) for mass ratios q = 1, 2, 3, 4 and 6. 

This procedure is not unique; we optimize a^*^ and Og ' 
over a relatively coarse grid and find reasonable results 
setting a(,{v)/v as a constant and 



a^{v) ^ (-5.828- 
ag(i^) = 184 i^. 



(36a) 
(36b) 



For this constant choice of aQ{v)/v^ we show in a subplot 
of Fig. 5 the a^(y)jv values constrained by the contours, 
as well as the quadratic fit of them. The error bars cor- 
respond to the width of the contours. More constraints 



imposed by new numerical waveforms may further lift the 
degeneracy between a^(y) and a^(y). We shall show in 
the next section that this model calibrated to the five nu- 
merical waveforms, even though not carefully optimized 
with a fine global grid, is good enough for detection with 
Advanced LIGO, and fair for parameter estimation pur- 
poses. 

Although a^(y)lv and a^iy^jv can not both be con- 
stants in this calibrated model, we shall emphasize that 
it does not imply that the physical 4PN and 5PN coef- 
ficients shall depend on v beyond the linear order. The 
optimal choice of a^(y) and a^{y) depends on other ele- 
ments of the dynamics, for instance, the Fade expression 
of A(r) or D{r), or the way we factorize or resum 
in the inspiral waveforms that enter the energy flux, etc. 
Therefore, it is possible that when a different EOB model 
is calibrated to the same set of numerical waveforms sub- 
ject to the constraint imposed by the self- force result, the 
optimal choice of a^{i>)/i> and aQ{v)/v are constants. It 
might be the case that a minor change in the dynamics 
of the EOB model being calibrated in this paper could 
bring all contours and the self-force constraint to cross 
at exactly the same point. 

Having calibrated the adjustable parameters of the 
EOB dynamics, we calibrate the adjustable parameters 
of the EOB waveform listed in Table I, notably the width 
of the comb At^^^^j^, a few higher-order FN terms in 
and 5ira (see Appendix B), and a pseudo QNM. For tech- 
nical reasons we do not include the adjustable parame- 
ters /921j P337 P44i and /955 in the energy flux (12), but 
only in the EOB modes (13). The energy flux being used 
in the dynamical evolution is therefore slightly different 
from the energy flux defined in Eq. (12). For the mass 
ratios and length of waveforms being considered in this 
paper, the fractional difference in the energy flux grows 
from ~ 10"'^ during low frequency to ~ 10~^ at merger. 
Such a difference in the dynamics generates a phase dif- 
ference at merger that increases from 0.02 rad to 0.35 
rad with mass ratio increasing from g = 1 to (7 = 6. In 
principle, to have an energy flux in the dynamics that is 
exactly consistent with the definition in Eq. (12), we need 
to calibrate those adjustable parameters in p^m as EOB- 
dynamics parameters, together with 05 and og. Calibrat- 
ing a large number of EOB-dynamics parameters has two 
consequences: the calibration becomes computationally 
expensive and all parameters become highly degenerate. 
Because of these technical difficulties and the relatively 
small fractional difference in the flux, we choose not to 
include the adjustable parameters in pim in the energy 
flux. 

We determine the width of the comb Atf^^^j^ by re- 
quiring the best agreement between EOB and numerical 
/i£„i around merger and ringdown. We find that Atf™^^jj 
depends moderately on the mass ratio, and we can as- 
sume At^™^j,jj being a constant. Specifically, we obtain 
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FIG. 6: For the equal-mass case, we compare the numerical- 
relativity and calibrated EOB (2, 2) mode. The top panels 
show the real part of numerical and EOB /122, the bottom 
panels show amplitude and phase differences between them. 
The horizontal axis is the retarded time in the numerical- 
relativity simulation. The left panels show retarded times 
t — r, = to 3850M, and the right panels show retarded 
times t - r, = 3850M to t - r* = 4070Af on a different 
vertical scale. The dotted curves are the numerical-relativity 
errors. 
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Calibrating the amplitude and phase of the EOB wave- 
form for the higher-order modes we find 
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As explained in Sec. JIB, since the iterative procedure 
that determines the NQC coefficients a^^^ usually takes 
4 to 5 steps to converge, we give fitting formulas of a^^^ 
as quadratic functions of the mass ratio to save compu- 
tational cost 



'■{v) ^ -4.559+ 18.76 1/- 24.23 
•(i^) = 37.68 - 201.51/ + 324.6 1^2^ 
'■{ly) = -39.6 + 228.9 1/- 387.2 1/^ 
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Finally, as discussed in Sec. II C, to improve the agree- 
ment of the (4, 4) and (5, 5) modes around merger, we 
introduce pseudo QNMs having 
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for all mass ratios considered in this paper. The fre- 
quency of these pseudo QNMs are about the frequency 
of the inspiral-plunge waveforms at the matching time 
and the decay time of these modes are about the 



match 



same as that of the first overtone of the physical QNMs. 



D. Comparing numerical and EOB (2, 2) mode 

In Figs. 6, 7 we compare the numerical-relativity and 
EOB (2,2) modes for mass ratios q = 1,2,3,4 and 6. 
We find that throughout the evolution the phase differ- 
ence is below ^ 0.1 rad. During the inspiral, the relative 
amplitude difference is within 2%, while during merger 
and ringdown it increases to within 12%. The numerical 
errors are also showed in the figures with dotted lines. 
We observe that during the inspiral the phase and am- 
plitude differences can be a factor of a few larger than 
the numerical-relativity error, but during the merger and 
ringdown they can be comparable or even smaller. As we 
shall see in Sec. IV, the mismatch between the numerical 
and EOB modes are consistently very small for detection 
purposes for Advanced LIGO, and the EOB modes are 
reasonably accurate for parameter-estimation purposes. 



E. Comparing numerical and EOB {I, m) (2, 2) 
modes 



In Figs. 8; 9 and 10, we compare the numerical and 
EOB subdominant modes ^133, /121 and /144 for the cases 
(7=1,3, 6. [For mass ratios (7 = 2,4 the plots look similar, 
so we do not show them.] During the inspiral, the numer- 
ical and EOB subdominant modes agree very well, simi- 
larly to the agreement we found for the /122 ■ This happens 
because the numerical frequencies ojim are well modeled 
by a simple multiple of the orbital frequency mfl. Dur- 
ing merger and ringdown, the agreement is very good for 
the ft.33 and h2i modes, i.e., comparable to the agreement 
of the /122 mode. Analogous performances hold for the 
other cases q = 2 and q = A. The numerical and EOB /144 
mode, however, show larger differences during ringdown. 
For instance, the phase difference increases to ^ 0.6 rad. 
There are two reasons for this less satisfactory result: (i) 
the larger errors in the numerical mode (4, 4) and (ii) 
the EOB QNM matching procedure that generates the 
ringdown part. Issue (i) spoils the numerical predictions 
of the fitting formulas for the (4, 4) mode (see Table III) 
which are essential to model the merger. Issue (ii) pre- 
vents modeling the ringdown phase of the ^144 with high 
accuracy (see Fig. 3 and discussions therein). Neverthe- 
less, since as seen in Figs. 1, the ft-44 amplitude is a few 
percent of the /i22 amplitude, the absolute error in ft.44 is 
generally smaller than the error with which we currently 
model the /i22 EOB mode. Therefore, the large difference 
between the numerical and EOB /144 is not the dominant 
source of systematic error in the gravitational polariza- 
tions. Since the mode comparison is very similar to 
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FIG. 7: Comparison between the (2, 2) modes of the numerical-relativity simulation and the calibrated EOB model for mass- 
ratios q — 2,3, 4, 6. In each sub-plot, the top panels show the real part of the /i22 mode and the bottom panels show the phase 
and amplitude differences between numerical and EOB waveform. The dotted curves are the numerical-relativity errors. 



that of the (4, 4) mode, except for an even larger phase 
difference of ~ 1 rad during ringdovifn, we do not show 
it for brevity. Modeling the /155 mode is difficult due to 
the same two issues discussed above, while on a more 
severe level. Nevertheless, we find in Sec. IV that there 
is substantial benefit in including this mode in the full 
polarization waveforms even though its modeling is not 
fully satisfactory. 

We point out that the special treatment of the (2,1) 
and (4,4) modes in Eq. (17), namely the replacement of 
D^"'''^'' with w^^^*^ ^V^f2 these modes, was necessary to 
improve the agreement of these modes with the numeri- 
cal waveforms. Eq. (17) was suggested by similar studies 
in the test- particle limit [60]. The reason is the follow- 
ing: as shown in Fig. 1 the amplitude of the numerical 
(2, 1) and (4, 4) modes reaches a peak quite after the peak 
of the (2,2) mode, i.e., they have large, positive Atp"^]^ 
values. If we want to impose the condition 1. listed in 



Sec. II B, the peak of the EOB mode should be moved 
to ipoak + ^^pTak- However, the EOB Newtonian ampli- 
tude is proportional to a power of the orbital frequency 
and the latter decreases to zero at the EOB horizon, thus 
the EOB amplitude drops to an extremely small value at 
^poak + ^*polk- By contrast, replacing vl = [rnnf with 
^/rn, we slow down the decay of these modes after t^^^^^, 
and we can successfully move the peak of the mode to 
*pcak + ^^pTak- Note that vl and 1/rn are exactly the 
same in the adiabatic Keplerian limit. This replacement 
introduces higher order non-adiabatic non-Keplerian cor- 
rections that have negligible effects on the inspiral wave- 
form and the energy flux. During plunge and merger, 
perturbative treatments break down and there is no a 
priori justification in describing the mode amplitude with 
a power of either w$ or l/r^. Equation (17) is adopted 
since combined with the NQC corrections, it is capable 
of reproducing the numerical-relativity waveforms dur- 
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FIG. 8: Comparison of the (4, 4) mode for mass-ratios q — 
1, 3, 6 between numerical and calibrated BOB model. Plotted 
data as in Fig. 7. 



FIG. 9: Comparison of the (2, 1) mode for mass-ratios g = 3 
(top panel) and q = 6 (bottom panel) between numerical and 
calibrated EOB model. Plotted data as in Fig. 7. 



ing merger. To make minimal adjustments to the New- 
tonian term of Eq. (16), we introduce this replacement 
only when needed, i.e. to the (2, 1) and (4,4) modes. 



IV. EFFECTUALNESS AND MEASUREMENT 
ACCURACY OF EOB WAVEFORMS 

In this section, we examine the effectualness and mea- 
surement accuracy of the EOB waveforms in matching 
the numerical-relativity waveforms. 

These investigations utilize the noise-weighted inner 
product between two waveforms hi, i = 1,2: 



{hl,h2) EE 4 5R 



Shif) 



df: 



(42) 



where hi{f) are the Fourier transforms of hi{t), and 
Sh{f) is the spectral density of noise in the detector. We 
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FIG. 10: Comparison of the (3, 3) mode for mass ratios g = 3 
(top panel) and q = 6 (bottom panel) between numerical and 
calibrated EOB model. Plotted data as in Fig. 7. 



choose one of the Advanced LIGO noise curves, named 
ZERO_DET_HIGH_P in Ref. [62]. 

Figure 11 shows the noise curve, and the amphtudes of 
the Fourier transforms of the numerical relativity wave- 
forms. For the binary's total masses considered, the NR 
waveforms start in band; to reduce artifacts from this, we 
taper the NR waveforms using the Planck-tapcr window 
function [63]. The width of the window function is set 
to the length of numerical relativity waveforms, which 
is about 0.5(Af/20MQ) seconds. The window function 
smoothly rises from to 1 in the first 0.0625 seconds 
and falls from 1 to in the last 0.0125 seconds. Further- 
more, whenever we compute quantities involving both a 
NR waveform and an EOB waveform (e.g., an overlap; 
see below), we restrict integration of the EOB waveform 
to the frequency for which numerical data is available. 

Figure 11 indicates the initial gravitational- wave fre- 
quency of the q — 1 and g = 4 simulations. The sim- 
ulation with q — 1 has the lowest initial gravitational- 




/(Hz) 



FIG. 11: Amplitude of the Fourier transform of the (2, 2) 
mode of the numerical-relativity waveforms, scaled to a to- 
tal mass M = 20Mq. We also plot the noise spectral den- 
sity of the Advanced LIGO detector. The two vertical lines 
mark the initial gravitational-wave frequency for the numeri- 
cal waveforms with q — 1 and q — 4 (lowest and highest start 
frequency of all considered waveforms). 



wave frequency, while the one with q = 4 has the largest 
initial gravitational-wave frequency. The numerical- 
relativity waveforms for mass ratios q ~ 1,2^ 3, 4, 6 have 
32,31,31,31, and 43 gravitational-wave cycles, respec- 
tively, from the initial frequency to the peak of the 
(2, 2) mode. Using the EOB model of this paper we 
estimate that for a total mass of 20Mq and mass ra- 
tios q = 1,2,3,4,6, there are 582,656,779,914,1184 
gravitational-wave cycles between 10 Hz and the start 
of the numerical-relativity simulations. These missing 
cycles, which decrease as the total mass of the binary 
increases, are not accounted for when computing mis- 
matches. 



A. Subdominant modes 

To investigate the importance of subdominant modes 
(Z,m) different from (2,2), we consider the gravitational 
waveform emitted from the binary into a given sky direc- 
tion {6, (p) (as measured relative to the orbital plane of 
the binary and note that is degenerate with the initial 
phase), given as 

h+{e, 0; t) - ihy, {9, -^Yi^rniO, (b) hem{t) . (43) 



Here -2Yim(0,(p) is the —2 spin- weighted spherical har- 
monic, and the summation on i and m is over the avail- 
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FIG. 12: The polarization waveform h+{6,(t);t) as emitted 
into sky direction 8 — tf) = n /3. Top panel: mass ratio 
g = 1; bottom panel: mass ratio q = 6. The sohd blue 
curve represents the numerical data, the red dashed curve 
the BOB model, and only late inspiral, merger and ringdown 
are shown. 



able NR or EOB modes. ^ In the remainder of this see- 
tion, we always ehoose 9 = 7r/3 and — tt/S, and assume 
a relative binary-detector configuration such that the de- 
tector is only sensitive to (i.e., an antenna pattern 
F+ = 1, Fx =0). A comprehensive study of arbitrary 
gravitational polarizations for all sky directions ^, is 
left to future work. 

Figure 12 shows the resulting waveforms 
/i4-(7r/3, 7r/3; for the mass-ratios q = \ and (7 = 6. 
This figure clearly shows that for (7 = 6, subdominant 
modes are more important, and one immediately expects 
that disregarding subdominant modes will have a larger 
effect for the 5 = 6 case. Let us now quantify these 
expectations. 



B. Effectualness 

The effectualness can be described by the mis- 
match (A^) between two time-domain waveforms h\ and 
/i2 (^0 7 '/'o I • Here, we consider all waveforms to be the 
-|- polarization evaluated in the direction 9 = (j) = it 
[see Eq. (43)]. We take hi to be the numerical relativity 
waveform at one of the simulated mass ratios for some 



given total mass M . The second waveform /12 is taken 
to be our calibrated EOB model, where we have explic- 
itly displayed the dependence of this waveform on some 
reference time and reference phase 0Oj as well as the 
masses mi and 7712 represented in the vector A of the 
parameters of the binary. 

The mismatch is given explicitly by [21] 



M = 1 



(/7.1,/t2(^O,0O, A)) 



where 



(44) 



(45) 



denotes the norm induced by Eq. (42). When search- 
ing for the signal waveform hi with the template 
h2{to, 00, A), the horizon distance is reduced by a factor 
Ai relative to searching with the perfect template hi, and 
the reduction in event rate is given by 1 — (1— A^)'^ sa iM.. 
Ideally, the maximization in Eq. (44) is over {to, </>0j ■^}; 
however, sometimes we choose to neglect maximization 
over A, as detailed below. 

Figure 13 presents several mismatch calculations for 
the equal-mass case. The solid lines compare the numer- 
ical relativity data to the leading (2, 2) mode of our EOB 
model. For these two curves, maximization of M is per- 
formed over {tQ,4>Q,X\. If the numerical waveform is rep- 
resented only by its (2, 2) mode, then the mismatches are 
very small, reaching ^ 10~^. However, if the 7-mode nu- 
merical waveform is used with all modes shown in Fig. 1, 
then the mismatch increases by about an order of mag- 
nitude, showing that subdominant modes are noticeable 
even for the 9 = 1 case. 

The two dashed curves in Fig. 13 use the calibrated 
EOB model with all five calibrated modes included. For 
these two curves, we maximize the mismatch only over 
{^o,</'o}, for technical convenience and to save compu- 
tational cost. Therefore the obtained mismatches are 
only upper bounds. We see that the 5-mode EOB model 
agrees significantly better with the 7-mode NR waveform 
than an EOB model utilizing only the (2, 2) mode. The 
line "5-mode NR vs. 5-mode EOB" compares NR with 
EOB when both waveforms contain only those five modes 
for which we calibrate the EOB model. 

In Figs. 13; 14 and 15 we have marked total mass 
IOOA/0 and 58.3Af0 for g = 1 and q = 6, respec- 
tively, the largest masses for stellar-mass binary black 
holes, assuming that the maximum black-hole mass is 
- 5OAf0.'^ More massive black holes are referred to 
as intermediate-ma.ss black holes (IMBHs). Their ex- 
istence and gravitational-wave event rates are more un- 
certain [67, 68]. 



In the numerical simulations, 7 modes (14 modes if we count 
m < 0) are extracted: (£,m) = (2,2), (2,1), (3,3), (3,2), (4,4), 
(5, 5), (6, 6). In the EOB model, 5 modes (10 modes if we count 
m < 0) modes are calibrated in Sees. HID and HIE: {£,m) = 
(2,2), (2,1), (3,3), (4,4), (5,5). 



As of today, the heaviest mass of a single black-hole source is 
23-34M0 [64, 65] , but considering the possibility of lower metal- 
licity we adopt as maximum mass of a black hole ~ 50Mq [66] . 
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Figure 14 presents the analogous calculations to Fig. 13 
for mass ratio q = 6. The mismatches are generally 
larger, owing to the more complex waveform of a (7 = 6 
binary. The numerical (2, 2) mode can still be fit by 
a (2,2) EOB waveform to 7W ~ 10~^. However, try- 
ing to represent the 7-mode NR waveform with only the 
EOB (2, 2) mode results in mismatches ^ 7% for to- 
tal mass M ^ 58.3M0 and above 10% for total mass 
M ~ 2OOM0. Thus, it is extremely important to ac- 
curately model higher-order modes when the merger and 
ringdown waveforms are in band and binary systems con- 
tain intermediate mass black holes with > 50Mq. So 
far, higher-order modes have been largely ignored in the 
analysis of real detector data (see for instance Ref. [33]), 
and by including a few dominant ones, the event rate 
or the horizon distance can be substantially increased, 
especially for high total masses. 

Inclusion of the higher-order modes reduces the mis- 
matches by a factor of 10 to ~ 5 x 10"'^ at low masses 
and ~ 0.01 at high masses, caused mainly by the error 
in modeling the ringdown waveforms of the higher order 
modes. The (3, 2) and (6, 6) modes that are not modeled 
or included in the EOB polarizations are not responsible 
for the comparatively large mismatch. To verify this, we 
show in these figures also the mismatches between the 
5-modc EOB polarizations and 5-mode numerical polar- 
izations that contain the same modes as the EOB ones.^ 
More diagnostic tests show that the error in modeling the 
(5, 5) mode is responsible for about half the increase in 
because the phase error of the (5,5) mode accumu- 
lates faster due to its higher frequency, and because it is 
larger than the phase error of other modes. 

As in Fig. 13, mismatches with the 5-mode EOB model 
are not optimized over A and represent therefore merely 
upper limits. 

We emphasize again that these mismatches mea- 
sure the difference between numerical-relativity and 
EOB waveforms only over the frequency band where 
numerical-relativity simulations are available (see 
Fig. 11). For large masses, say M = IOOMq, the initial 
frequency of the numerical-relativity (2, 2) mode for 
q = 6 is 15Hz, and the missing part between lOHz 
and 15Hz would likely modify the mismatches only 
marginally because of the steeply rising seismic noise 
wall toward low frequency. However, for M = 20Mq and 
(7 = 6, the initial frequency of the numerical-relativity 
(2, 2) mode is 70Hz which is in the frequency band of 
the detector. In this case we miss the portion of the 
numerical-relativity waveform between lOHz and 70Hz. 
The best way, though time-consuming, to address this 
gap problem is to produce longer numerical-relativity 
waveforms, or to conduct tests within the analytical 
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Sometime the mismatches increase by ~ 0.1% after we remove 
the (3, 2) and (6, 6) modes from the numerical-relativity polar- 
izations. This happens because the numerical precision of our 
code in optimizing over the initial phase is ~ 0.1% 
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FIG. 13: The mismatch versus binary total mass for q — I us- 
ing Advanced LIGO noise curve. Mismatches with the (2, 2) 
mode of the EOB waveform are optimized over {to,(l>o,X} 
whereas mismatches with the 5-mode EOB model are opti- 
mized over {to,(j}o} only and represent an upper bound. The 
vertical line represents the maximum total mass for stellar- 
mass black-hole binaries, assuming a maximum black-hole 
mass of 50A/q . The range of total masses on the right of the 
vertical line refer to intermediate-mass black-hole binaries. 



models to assess their reliability below a certain fre- 
quency [21, 59, 69-72]. With the caveat of this frequency 
gap, we conclude that our calibrated EOB model is 
sufficient for detection purposes. 



C. Measurement accuracy 

We are now interested in the question of whether the 
EOB polarizations are accurate enough to be used in data 
analysis for measurement purpose. We adopt as accuracy 
requirement for measurement the one proposed by [73, 
74] 



\\Sh\\ < e. 



(46) 



where Sh = h^'~'^{t) — h^^{t) is the error in modeling 
the numerical waveforms, and e < 1 incorporates a safety 
factor [71] or the effect of a detector network [59]. The 
left hand side in Eq. (46) increases proportionally with 
the signal-to- noise (SNR) ratio and we calculate 

the upper bound of the SNR"^ = SNR/e that satisfies 
Eq. (46). For any SNR below this upper bound, the 
EOB waveforms or polarizations are accurate enough for 
measurement purposes, i.e., accurate enough not to gen- 
erate any systematic bias that is larger than statistical 
errors in estimating the physical binary parameters. The 
upper bound SNR"*^ is therefore a very strict accuracy re- 
quirement. Unlike the well known fact that good phase 
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FIG. 14: Mismatch calculation for mass ratio q = 6. All 
details as in Fig. 13. This figure and Fig. 13 use the same 
vertical axis to ease comparisons between them. 



agreement is critical for obtaining good effectualness, to 
get high (upper bounds on) SNR°^, both the amplitude 
and the phase of the templates must agree very well with 
those of the exact waveforms. 

In Fig. 15, for mass ratios q = I and q = 6, we show 
the upper bound of the SNR'^^ as a function of the total 
mass. These SNR°*^ are calculated for a single Advanced 
LIGO detector. In the q = 1 case, the 5-mode EOB 
polarizations match the 7-mode numerical polarizations 
accurately enough for any SNR''^ < 24 when the total 
mass IS below lOOMg, and for any SNR'=^ < 19 when 
the total mass is below 2OOM0. These upper bounds of 
the SNR'^* may not seem impressive at the first sight. 
However, there is a significant improvement if we com- 
pare them with the upper bounds obtained when only 
the EOB (2,2) mode is used which are also shown in 
Fig. 15. Note that the right vertical axis in Fig. 15 shows 
||(5/i||/||/i||. Once the EOB and NR waveforms are aligned 
at low frequency, we do not allow further time or phase 
shift in calculating Sh. This differs from the practice of 
Refs. [59, 71] which plot minimized over time- 

and phase-shifts. 

In the q = 6 case, the upper-bound SNRs are lower, 
because of the larger error in modeling the higher-order 
modes and the relatively high contribution to the SNR 
from higher-order modes for such an asymmetric binary. 
For stellar-mass black holes, M < 58. SM©, the 5-mode 
EOB polarizations are accur ate for SNR'^* < 10, and 
when the total mass is below 200 AIq, the 5-mode EOB 
polarizations are accurate for SNR'^* < 5. The higher 
order modes in the EOB model, especially their ringdown 
part, clearly needs better modeling. Nevertheless, the 
improvement from using only the EOB (2,2) mode as 
templates is significant. 



FIG. 15: The upper bound SNR'=" = SNRe from the mea- 
surement accuracy requirement (46) versus binary total mass 
for q = 1 and 6 using Advanced LIGO noise curve. For any 
SNR''^ below the curves, the EOB polarizations are accurate 
enough to avoid a systematic bias that is larger than statistical 
errors when estimating the binary parameters. The horizon- 
tal line indicates the single detector SNR for Advanced LIGO 
which is 8. The two vertical lines represent the maximum 
total mass for stellar-mass binary black holes with maximum 
black-hole mass of 50Mq and q = 1 and q = 6. The range 
of total masses on the right of those vertical lines refer to 
intermediate-mass black-hole binaries. The vertical axis on 
the right shows 



In closing, we mention a few important caveats of our 
analysis: First, all norms and mismatches computed for 
Figs. 13-15 are performed only over those frequencies for 
which numerical data is available. At lower total mass, 
more gravitational- wave cycles lie in the LIGO frequency 
band, so this restriction becomes more severe. Therefore, 
at low masses, our analysis might yield an increasingly 
overoptimistic view. 

On the other hand, the requirement (46) on the ac- 
curacy measurement may be unnecessarily stringent. In 
fact, although Fig. 15 would say that for certain mass 
ratios and total masses, the EOB waveforms will intro- 
duce biases in the binary parameters, these biases may 
affect binary parameters that have little astrophysical 
relevance — for example the gravitational-wave phase at 
coalescence. 

Finally, we emphasize again that the results in this 
section are just a first step into the problem of including 
higher-order modes, since we investigated only one geo- 
metrical configuration. A comprehensive investigation of 
the gain in effectualness and measurement accuracy, as 
well as the impact on estimating the source parameters, 
when including these rather accurate higher-order modes 
will be presented in a separate paper. 
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V. RELATION TO PREVIOUS WORK 

The EOB model we consider here differs from the non- 
spinning EOB model employed by Buonanno et al. [23] 
in the handling of the radiation-reaction sector. Refer- 
ence [23] adopted a Pade-resummed radiation-reaction 
force and energy flux [22, 75], while here we adopt the 
factorized- waveform energy flux of Refs. [26, 30, 76]. Ref- 
erences [25, 26] found that when generalizing the Pade- 
resummed flux to the spin case the agreement with the 
numerical energy flux is not very satisfactory. Thus, 
we concluded that the nonspinning EOB model with 
Pade-resummed flux is not a very good candidate for the 
generic spin EOB model. 

By using the SpEC numerical merger (2, 2) mode with 
mass ratio g = 1, and inspiral (2, 2) mode with q = 2,3, 
Buonanno et al. [23] calibrated the 4PN parameter in the 
EOB radial potential A(r) and the parameter iipoie in the 
Padc-rcsummcd energy flux, such that the EOB model 
could be also used outside the range of binary masses 
employed to calibrate it. We find that when comparing 
to the SpEC merger (2, 2) mode of this paper with mass 
ratios q = 1,2,3, the EOB (2,2) mode of Ref. [23] have 
maximum phase difference until merger of 0.12, 0.22, and 
0.09 rad, respectively. In contrast, the model calibrated 
here results in smaller or comparable phase differences of 
0.04, 0.08 and 0.12 rad, respectively. Numerical data for 
g = 4, 6 were not available for the calibration in Ref. [23]. 
Comparing the new numerical data available now with 
the model from Ref. [23], we find phase differences of 0.43 
and 1.8 rad, respectively, for g = 4 and 6 (in contrast, our 
new model results in 0.12 and 0.15 rad, respectively; see 
Fig. 7). Those phase differences are obtained by using 
the low- frequency alignment procedure of Eq. (34). If 
wc were adopting the two-pinching frequency procedure 
of Ref. [31], we would obtain 0.065 rad with pinching 
frequencies MD. = 0.052 and = 0.3 for q ^ 4, and 0.18 
rad with pinching frequencies Mfl = 0.056 and = 0.15 
for q = 6. 

Before calibration, the radiation-reaction sector of the 
EOB model used here almost coincides^ with the non- 
spinning EOB model used in Damour and Nagar [31] in 
the radiation-reaction sector, but only in its uncalibrated 
version. In fact, the adjustable parameters used here dif- 
fer from the ones used in Ref. [31]. Moreover, in this pa- 



Whereas we use the factorized waveforms with pi^ny Ref. [31] 
employs factorized waveforms where (p22)^ is traded with its 
Taylor expanded form and then Pade resummed [30]. In the 
factorized waveforms, Ref. [31] also adopts a different value for 
the constant in the tail term ro = 2, a different odd-parity source 
term S^g (r, pr, , p$ ) = p<j /vq , and a different Ni,„ given in Eq. 
(5) of Ref. [31]. Furthermore, while we include radiation reaction 
force in Eq. (10c), i.e., in the equation of motion of p^, , Ref. [31] 
does not. Finally, Ref. [31] trades the Keplerian velocity with 
the non-Keplerian velocity in T^mi Plm but not Stm, and does 
not include the higher-order PN terms computed in Ref. [26] . 



per we also introduce adjustable parameters in the phase 
of some of the {i, m) modes, and in some cases also in the 
factorized amplitude. Furthermore, we also modify the 
leading Newtonian term for the (2,1) and (4,4) modes 
[see Eq. (17)]. 

By using the SpEC numerical-relativity merger (2, 2) 
mode with mass ratio q = I, and merger (2,2) mode 
with q = 2,4 from the Jena numerical-relativity group, 
Damour and Nagar [31] calibrated the 4PN and 5PN pa- 
rameters in the radial potential A{r). Wc find that when 
comparing to the SpEC (2, 2) modes of this paper with 
mass ratios q = 1,2, 4, the EOB (2, 2) mode of Ref. [31] 
has maximum phase difference until merger of 0.25 rad, 
0.36 rad, 1.32 rad, respectively, when aligning at low fre- 
quency and 0.05 rad, 0.11 rad, 0.25 rad when using the 
two-frequency pinching procedure. In the case of mass 
ratios q = 3 and 6 where numerical waveforms were not 
available for calibration, the phase differences increase to 
0.93 rad and 2.3 rad, respectively, when aligning at low 
frequency and 0.17 rad, 0.65 rad, when using the two- 
frequency pinching procedure. 

We notice that once we have calibrated the EOB model 
to a set of numerical-relativity waveforms using the low- 
frequency alignment procedure of Eq. (34), the phase dif- 
ference with numerical-relativity waveforms is not sensi- 
tive to the alignment method — for example we find the 
same phase differences when using the two-pinching fre- 
quency procedure. By contrast, EOB waveforms cali- 
brated with the two-frequency pinching procedure do not 
have the same phase difference with numerical-relativity 
waveforms when aligning them at low frequency. We de- 
duce then that the low-frequency alignment procedure is 
more robust. 

Recently, the LIGO/ Virgo detectors have completed 
the first search for gravitational waves from high-mass 
black holes using analytical inspiral, merger and ring- 
down templates [33]. They employed the EOB model 
calibrated to numerical waveforms from NASA-Goddard 
of Ref. [20]. The mismatches between the EOB and 
NASA-Goddard waveforms computed in Ref. [20] are less 
than 3% when the total mass is 25-99M0, which is the 
mass range used in the LIGO/Virgo search [33]. These 
mismatches were derived maximizing only on the time 
and phase, but not on the binary parameters. We find 
that when comparing the EOB (2,2) mode of Ref. [20] 
with the EOB model developed in this paper, the phase 
and amplitude differences at merger arc at most of 5 
rad and 20%, respectively, when q ~ 1,2,3,4, and 5.8 
rad for q = 6. The EOB higher-order modes (2,1), 
(3,3), and (4,4) were also calibrated for the first time 
in Ref. [20] using the NASA-Goddard numerical wave- 
forms. In this case the matching between the inspiral- 
plunge and merger-ringdown waveforms was chosen at 
the same point in time for all the modes. Those higher- 
order modes were not employed in the search of Ref. [33] . 
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VI. CONCLUSIONS 

The first search for gravitational waves from nonspin- 
ning high-mass binary black holes (M = 25-99Mq) with 
LIGO and Virgo detectors has been recently completed, 
setting astrophysically meaningful upper limits [33] . The 
search has used for the first time templates which include 
inspiral, merger and ringdown. Those templates were 
built by combining numerical-relativity and analytical- 
relativity results, cither through the EOB waveforms 
of Rcf. [20] (see also the most recent improvements in 
Refs. [23, 25, 28, 29, 31, 32]) or the phenomenological 
merger-ringdown waveforms proposed in Refs. [77] (see 
also Ref. [78] for an updated version). 

In this paper we have built on Refs. [19, 20, 23, 25, 28- 
31], and have improved further the EOB model taking 
advantage of highly accurate numerical-relativity simu- 
lations with mass ratios q = 1,2, 3, 4, 6 from the Caltech- 
Cornell-CITA cohaboration [41]. 

By extracting several numerical-relativity quantities, 
such as the mode's amplitude and its second time deriva- 
tive at the peak, as well as the frequency and its first time 
derivative at the peak, we have improved the agreement 
of numerical and EOB phase and amplitude very close to 
merger, reproducing important features of the numerical 
simulations — for example the fact that whereas the (2, 2) 
mode peaks at the same position of the EOB light ring, 
the higher-order modes peak at late times [20, 32, 51- 
53, 60]. 

We have found that the (2, 2) mode can be calibrated 
with M < 0.1% for M = 20-200 A/q and mass ra- 
tios q = 1,2,3,4,6, using only three EOB adjustable 
parameters, notably the 4PN and 5PN terms, as and 
ae, in the radial EOB potential, and the width of the 
comb, ^t'^i.tch (^66 Table I). We have also found that 
the strongest subdominant modes (2, 1), (3,3), (4,4) and 
(5, 5), can be successfully calibrated by including for each 
mode four EOB adjustable parameters, specifically the 
3PN terms in p2i,P33, P44, P55, the 2.5PN or 3.5PN 
terms in (S2ii'533, and S^^, the width of the comb 

A^match: Adtch: A^match- ^*match' ^ud, in SOmC CaSCS, a 

pseudo QNM (see Table I). The reason of introducing 
more parameters for higher modes rests on the fact that 
those modes are known at PN orders lower than the (2, 2) 
mode. Furthermore, to achieve this very good agreement 
of the modes' phase and amplitude we have also used the 
information from numerical-relativity about the peak's 
amplitude and frequency of Table III. and the final 
masses and spins of Eq. (29). These data determine the 
complex amplitudes entering the merger-ringdown wave- 
form (28), and the NQC coefficients in Eq. (22). 

When investigating the cffectualness for detection pur- 
poses, we have found that the numerical-relativity po- 
larizations containing the strongest seven modes have 
a maximum mismatch of 7% for stellar-mass binaries, 
and 10% for intermediate mass binaries, when only the 
EOB (2,2) mode is included for q = 1,2,3,4,6 and bi- 
nary total masses 20-200 Hz. However, the mismatches 



decrease when all the four subleading EOB modes cali- 
brated in this paper are taken into account reaching an 
upper bound of 0.5% for stellar-mass binaries, and 0.8% 
for intermediate mass binaries. Event rates or horizon 
distances can be substantially increased, especially for 
high total masses, if those subleading modes were in- 
cluded in gravitational- wave searches [33] . 

We have also emphasized that when computing the 
mismatches, we do not attach any PN waveforms to the 
numerical-relativity waveforms, because we do not want 
to introduce any error due to the procedure of build- 
ing hybrid PN-numerical waveforms. Moreover, for bi- 
naries with low total mass — say 20-100Mq — many more 
gravitational-wave cycles than the ones of the numerical 
simulations used in this paper are in band. Our mis- 
matches do not take into account these missing cycles. 
As a consequence if the EOB model were used to search 
for signals of length larger than the one of the numerical 
waveforms employed here, the mismatches could become 
larger. 

We have also studied the measurement accuracy of the 
EOB model using the accuracy requirement proposed in 
Ref. [74]. Using one single Advanced LIGO detector, we 
have determined the SNRs below which the EOB polar- 
izations are accurate enough that systematic biases are 
smaller than statistical errors. Unlike the well known 
fact that good phase agreement is sufficient for obtain- 
ing good effectualness, to get high upper-bound SNRs, 
both the amplitude and the phase of the templates must 
agree very well with those of the exact waveforms. Since 
higher-order modes have non-negligible contribution for 
large mass ratios, and those modes have the largest am- 
plitude errors, we have found that the upper-bound SNRs 
are lower for the most asymmetric systems. We stress 
again the relevance of modeling the higher-order modes, 
because using only the (2, 2) mode would decrease signif- 
icantly the upper-bound SNRs. However, it is worth to 
note that the accuracy requirement that we used may be 
too stringent since by itself it does not say which of the 
binary parameters is going to have biases [79]. It could 
turn out that the biased parameters have little relevance 
in astrophysics or tests of general relativity. 

Finally, we have used rather long numerical-relativity 
waveforms, in particular the case q = 6 has forty 
gravitational- wave cycles before merger. Confirming pre- 
vious studies [21, 59, 69-72] we have found, that espe- 
cially for large mass ratios, the addition of more cycles 
at low frequency does affect the accuracy of the EOB 
waveforms (as any other PN waveforms) which were cal- 
ibrated to a shorter number of cycles. So, we had to re- 
calibrate the EOB adjustable parameters to achieve very 
small phase errors around merger. This re-calibration 
is crucial for parameter estimation, but not for detec- 
tion, and we expect to do it again in the future when 
longer or more accurate numerical-relativity waveforms 
for asymmetric systems will become available. More- 
over, we plan to improve the matching procedure from 
inspiral-plunge to merger-ringdown, since the majority 
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of the phase and amplitude error is accumulated dur- 
ing this transition, especially for higher-order modes. Of 
course, further improvement of the EOB higher-order 
modes also depends on the availability of sufficiently ac- 
curate numerical-relativity data especially during the last 
stages of inspiral, merger and ringdown. 
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FIG. 16: Investigation of the (3, 2) mode during ringdown 
for q = 1 and 6. The top panels show the mode itself, 
and the lower panels split the mode into amplitude and 
frequency. The continuous lines represent the numerical- 
relativity modes. The dotted lines are a fit to eight QNM 
modes of {I = 3, m = 2,n = 0, 1, . . . , 7). The dashed lines 
are a fit to eight QNMs of = 3, m = 2, n = 0, 1, . . . , 4) and 
= 2, m = 2, n = 0, 1, 2). The horizontal axis is the retarded 
time in the numerical-relativity simulation. 



Appendix A: The {i, m) = (3, 2) mode 

Reference [19] found that the numerical-relativity (3, 2) 
mode contains QNMs with (£, to) = (3, 2) and (£, to) = 
(2, 2). This beating of QNMs with the same m but differ- 
ent i arises in the transformation from the spin- weighted 
spheroidal harmonics (which are eigenmodes of the ra- 
diation generated by the perturbed final black hole) to 
the spin-weighted spherical harmonics that are used to 
decompose the multipolar waveform. This is a general 
feature of modes with i > 2 and m < £, and since the 
(3, 2) mode is the dominant one among such modes, we 
discuss its modeling and possible calibration in this sec- 
tion. 

First, we confirm the result of Ref. [19] that the 
ringdown portion of the (3, 2) mode can be accurately 
modeled by a linear superposition of QNMs with both 
{e,m) = (3,2) and {i,m) = (2,2) modes. The result 
of Ref. [19] was restricted to the equal-mass case. Here 
we extend this analysis to unequal-masses and add more 
overtones. We fit the numerical ringdown mode with ei- 
ther (i) a set of eight QNMs with (^ 3, to = 2, n 
0, 1, . . . , 7) or (ii) a set of eight QNMs with {i = 3, to = 
2,n = 0,1,...,4) and (£ = 3,TO = 2,n = 0,1,2). For 
mass ratios q = 1 and 6, we compare the fitting results 
with numerical waveforms in Fig. 16. The results for 
other mass ratios are similar. The very different per- 
formance of sets (i) and (ii) shows clearly that the nu- 
merical ringdown mode does contain contributions from 
{i,m) = (2,2) QNMs. 
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FIG. 17: Amplitude and frequency comparisons between full 
"NR" (3, 2) mode and "NR-f QNM" (3, 2) mode generated by 
attaching the QNMs to the inspiral-plunge numerical wave- 
form. The "NR-I-QNM" waveform is generated by attaching 
the following set of eight QNMs (^ = 3, m = 2, n = 0, 1, . . . , 4) 
and (^ = 2, m = 2, n = 0, 1, 2). Note that the amplitudes have 
been multiplied by a factor of 20, so that they are more visi- 
ble. The horizontal axis is the retarded time in the numerical- 
relativity simulation. 
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Second, for the case q = 6, we explore in Fig. 17 
the possibility of modeling the numerical ringdown (3, 2) 
mode as a linear superposition of the QNMs of set (ii), 
using the EOB matching procedure of Sec. II C. We ob- 
serve that although the beating of the {£, m) = (3, 2) and 
{£, m) — (2, 2) QNMs reproduce qualitatively well the os- 
cillations in the ringdown amplitude and frequency, the 
behavior of the two curves "NR" and "NR-I-QNM" is very 
different. Therefore, we deduce that the current match- 
ing procedure which use the information of a small seg- 
ment of inspiral-plunge waveform before the peak of the 
amplitude, does not provide us with the correct QNMs 
coefficients when different i values are present. Since 
this matching procedure is applied to build the calibrated 
EOB modes, we find that the same problem affects the 
calibration of the EOB (3, 2) mode. We postpone to the 
future the solution to this important issue. 



Appendix B: Quantities entering the EOB 
gravitational modes 
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Following Ref. [30], we introduce the velocity param- 
eter, V = {^l H'"^'^^Y/'^ . The explicit expressions of the 
function Sim then read [26, 30] 
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Note that the 2.5PN and 3.5PN coefficients (5^i\4? and 
(7) 

(544 ill Eqs. (Bib), (B2a) and (B3a) are not known in 
PN theory. They are determined by calibrating the EOB 
to numerical-relativity waveforms. Their explicit expres- 
sions for these parameters are given by Eqs. (39). 

The following quantities enter the Newtonian modes in 
Eq. (16) 
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The odd-parity modes p^^^ and even-parity modes read [26, 30] 
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7187914 
15 526 875 



-a 1 



8 !/2 _ 626 1/ + 319 2 
390 (2 1/ - 1) " 

273 Z/3-861I/2 



P65 



P64 



P63 



P62 = 



L 

Pel 



31877 4 
304200 

- 602 1/ - 106 



84(5t/2-5t/+l) " 
220z/3-910zy2_|_ 838j,_ig5 ^ 

144 (3 1/2 -4 J/ +1) ^ 
133j/3-581i/2 + 462t/-86 , 



84(5i/2 - 5j/+ 1) 
156 - 750 i/2 + 742 i/ - 169 

144(3!/2_4j/+1) 
49 1/3 - 413 1/2 + 378 1/ - 74 2 

84(5^/2- 5;/+ 1) 
124 1/3 - 670 1/2 + 694 1/ - 161 

144 (3t/2-4j/+l) 



1025 435 
659 736 



476 887 
659 736 



817991 
3 298 680 



P77 = 1 



P76 = l 



P75 



P74 = l 



P73 



L 
P72 



P71 = 1 



13801/3 - 49631/2 + 4246J/ - 906 2 

714 (3i/2 - 4i/ + 1) 
61041/4 - 29351i/3 + 378281/2 - 16185;/ + 2144 . 



804^3 



1666 (7z/3 - 14i/2 + 71/ - 1) 
3523t/2 + 3382;/ - 762 , 



714 (3^2 _ _^ j^) fi' 
410761/4 - 217959^3 ^ 298872z/2 - 131805j/ + 17756 

14994 (7j/3 - 14z/2 + 7^-1) 
420J/3 - 25631/2 + 2806:/ - 666 2 

714 (3^2 _4j,_^ 1) > 

32760J/4 - 1902391/3 + 273924^2 „ i23489i/ + 16832 



228i/3 



14994 (7j/3 - 14^2 + _ 1) 
2083^2 _,_ 25181/ - 618 , 



714 (3^2_4j,^j^) 



(B12a) 
(B12b) 
(B12c) 
(B12d) 
(B12e) 

(B13a) 
(B13b) 
(B13c) 
(B13d) 
(B13e) 
(B13f) 

(B14a) 
(B14b) 
(B14c) 
(B14d) 
(B14e) 
(B14f) 
(B14g) 



Here we use eulerlog,„(w^) 



IE 



log 2 + logm 



are determined by calibrating the EOB to numerical- 



1/2 log Vq, with 'jE being the Euler constant. Note that relativity waveforms. The calibrated expressions are 
the 3PN coefficients p^^^ and pf^ in Eqs. (B9b), given by Eq. (38). 
(BlOa) and (Blla) are not known in PN theory. They 
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